function X = rk4(f, tspan, X0)
% rk4.m ——— 简单的四阶 Runge–Kutta 积分器
% 用法： X = rk4(f, tspan, X0)
%   f     — ODE 函数句柄，格式 @(t,X)
%   tspan — 时间向量，例如 0:dt:T
%   X0    — 初始状态列向量

    n = numel(tspan);
    dim = length(X0);
    X = zeros(dim, n);
    X(:,1) = X0;
    for i = 1:n-1
        h = tspan(i+1) - tspan(i);
        t = tspan(i);
        k1 = f(t,           X(:,i));
        k2 = f(t + h/2,     X(:,i) + h/2 * k1);
        k3 = f(t + h/2,     X(:,i) + h/2 * k2);
        k4 = f(t + h,       X(:,i) + h   * k3);
        X(:,i+1) = X(:,i) + h/6 * (k1 + 2*k2 + 2*k3 + k4);
    end
end
